;+
; NAME:
;       OMNI_ASTROCUT_FIT
;
; PURPOSE:
;       To help me graduate.
;
; CATEGORY:
;       AAAAAAAAAAAAAAAHHHHHHHHHHHHHH!!!!!!!!!!!!
;
; CALLING SEQUENCE:
;       HHH
;
; INPUTS:
;       NONE
;
; OPTIONAL INPUTS:
;       NONE
;
; KEYWORD PARAMETERS:
;       NONE
;
; OUTPUTS:
;       NONE
;
; OPTIONAL OUTPUTS:
;       NONE
;
; COMMON BLOCKS:
;       NONE
;
; MODIFICATION HISTORY:
;
;       Created:  03/13/14, TPEB -- Initial version.
;
;-

PRO OMNI_ASTROCUT_FIT
  
  COMPILE_OPT IDL2, LOGICAL_PREDICATE
  
  COMMON OMNI_CONFIG, conf, mw, local, dpdfs, ancil, fmt, conffile
  COMMON MASS_BLOCK, pvec, s, iflux, ieflux, m0, hnuk, kappa, nkp
  
  ;; List of everything the plotting subroutine needs
  COMMON CUT_PLOT, yr, dyr, mmin, mmax, alpha, mconf, ndloop, sunsym
  
  ;; Timing for length of run
  t_start = systime(1)
  nreal   = 1
  ndloop  = nreal
  
  ;; Read in survey-info, galactic-params, & dpdf-params config files
  conf     = omni_load_conf(cfile)
  mw       = omni_read_conffile('./conffiles/galactic_params.conf')
  dpdfs    = omni_read_conffile('./conffiles/dpdf_params.conf')
  mconf    = omni_read_conffile('./conffiles/mass_deriv.conf')
  d        = dindgen(dpdfs.nbins)*dpdfs.binsize + dpdfs.binstart
  sampletd = n_elements(sampletd) ? sampletd GE 1 : mconf.sampletd
  
  status = OMNI_MASS_SETUP(CONFFILE=cfile, VERBOSE=verbose, SILENT=silent, $
                           TYPSTR='astrocut')
  IF ~status THEN RETURN
  
  ;; Read in the SURVEY catalog 
  restore,'./local/'+conf.survey+'_map_locations.sav',VERBOSE=verbose
  
  ;; Limit SURVEY structure to the items saved in the above line
  surv_ind = WHERE_ARRAY(s.cnum, survey.cnum, n_survey)
  IF n_survey NE n_elements(s) THEN message,'Rum now, Think later!!!'
  survey = survey[surv_ind]
  restore,'./local/'+conf.survey+'_velocities.sav',/ver
  v = v[surv_ind]
  
  ;; Restore the PHYSICS variable
  restore,'./masses/save_files/physics.sav',/ver ; Restore the physics
  undefine,mmin,mmax,rmin,rmax,nmin,nmax         ; Clear these out; don't need
  
  ;; Set maximum mass to use in computation of mfns
  m_max = mconf.mmax
  message,'The maximum mass to be used for fitting functional forms is '+$
          string(m_max,format="(I0)")+' M_sun.',/inf
  print,''
  
  
  ;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  ;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  ;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  ;; Now, we are ready to conquor the world!  Hail Ceasar!
  
  ;;******************************************************************
  ;; Heliocentric Distance Cuts
  print,''
  print,'--------------------------------'
  print,'Heliocentric Distance'
    
  myps,'./masses/plots/astrocut_dsun.eps',xsize=10
  print,!x.omargin,!y.omargin
  !x.omargin = [-7,1]
  !y.omargin = [-1,-3]
  multiplot,[2,1],xgap=0.002
  sunsym = cgSymbol('sun')
  
  i1 = where(pvec.stat.duse[0]/1.d3 LE 6.5, n1)
  i2 = where(pvec.stat.duse[0]/1.d3 GT 6.5, n2)
  nt = n1+n2
  print,nt
  
  ;; Compute PL fit
  mass = 10.^phys[i1].mass[0]
  mmax = m_max
  alpha = OMNI_MFN_PL(mass, mmax, mmin)
  
  OMNI_PLOT_CUTS, mass, /first
  
  al_legend,/top,/left,box=0,textcolor='red',$
            [string(cgSymbol('alpha'),alpha,format="(A0,' = ',F0.2)"),$
             'Nearby Objects','d'+sunsym+' < 6.5 kpc']
  
  multiplot                     ;,/doyaxis
  ;; Second Panel *****************************
  
  ;; Compute PL fit
  mass = 10.^phys[i2].mass[0]
  mmax = m_max
  alpha = OMNI_MFN_PL(mass, mmax, mmin)
  
  OMNI_PLOT_CUTS, mass, /last
  
  al_legend,/top,/left,box=0,textcolor='red',$
            [string(cgSymbol('alpha'),alpha,format="(A0,' = ',F0.2)"),$
             'Distant Objects','d'+sunsym+' > 6.5 kpc']
  
  myps,/done,/mp
  
  
  ;;******************************************************************
  ;; Galactocentric Radius Cuts
  print,''
  print,'--------------------------------'
  print,'Galactocentric Radius'
    
  omni_lbd2rz,pvec.glon,pvec.glat,pvec.stat.duse[0],Rgal,zgal,$
              XGAL=xgal,YGAL=ygal,THETA=theta
  
  myps,'./masses/plots/astrocut_rgal.eps',xsize=10
  print,!x.omargin,!y.omargin
  !x.omargin = [-7,1]
  !y.omargin = [-1,-3]
  multiplot,[3,1],xgap=0.002
  sunsym = cgSymbol('sun')
  
  i1 = where(Rgal/1.d3 LE 5.5, n1)
  i2 = where(Rgal/1.d3 GT 5.5 AND Rgal/1.d3 LE 8.5, n2)
  i3 = where(Rgal/1.d3 GT 8.5, n3)
  nt = n1+n2+n3
  print,nt
  
  ;; Compute PL fit
  mass = 10.^phys[i1].mass[0]
  mmax = m_max
  alpha = OMNI_MFN_PL(mass, mmax, mmin)
  
  OMNI_PLOT_CUTS, mass, /first
  
  al_legend,/top,/left,box=0,textcolor='red',$
            [string(cgSymbol('alpha'),alpha,format="(A0,' = ',F0.2)"),$
             'Ring / Sct-Cen Arm','R!dgal!n < 5.5 kpc']
  
  multiplot                     ;,/doyaxis
  ;; Second Panel *****************************
  
  ;; Compute PL fit
  mass = 10.^phys[i2].mass[0]
  mmax = m_max
  alpha = OMNI_MFN_PL(mass, mmax, mmin)
  
  OMNI_PLOT_CUTS, mass
  
  al_legend,/top,/left,box=0,textcolor='red',$
            [string(cgSymbol('alpha'),alpha,format="(A0,' = ',F0.2)"),$
             'Sgr / Local Arm','5.5 kpc < R!dgal!n < 6.5 kpc']
  
  multiplot                     ;,/doyaxis
  ;; Third Panel *****************************
  
  ;; Compute PL fit
  mass = 10.^phys[i3].mass[0]
  mmax = m_max
  alpha = OMNI_MFN_PL(mass, mmax, mmin)
  
  OMNI_PLOT_CUTS, mass, /last
  
  al_legend,/top,/left,box=0,textcolor='red',$
            [string(cgSymbol('alpha'),alpha,format="(A0,' = ',F0.2)"),$
             'Per / Outer Arm','R!dgal!n > 8.5 kpc']
  
  
  myps,/done,/mp
  
  ;;******************************************************************
  ;; Physical Radius Cuts
  print,''
  print,'--------------------------------'
  print,'Physical Radius'

  ires = where(finite(s.rad),nres)
  iunr = where(~finite(s.rad),nunr)
  prad = phys[ires].prad[0]
  
  myps,'./masses/plots/astrocut_prad.eps',xsize=10
  print,!x.omargin,!y.omargin
  !x.omargin = [-7,1]
  !y.omargin = [-1,-3]
  multiplot,[2,1],xgap=0.002
  sunsym = cgSymbol('sun')

  i1 = where(prad LE alog10(0.125), n1)
  i2 = where(prad GT alog10(0.125) AND prad LE alog10(1.25), n2)
  i3 = where(prad GT alog10(1.25), n3)
  nt = n1+n2+n3+nunr
  print,nt
  
  ;; ;; Second Panel *****************************
  
  ;; Compute PL fit
  mass = 10.^phys[i2].mass[0]
  mmax = m_max
  alpha = OMNI_MFN_PL(mass, mmax, mmin)
  
  OMNI_PLOT_CUTS, mass, /first
  
  al_legend,/top,/left,box=0,textcolor='red',$
            [string(cgSymbol('alpha'),alpha,format="(A0,' = ',F0.2)"),$
             'Clump','0.125 pc < R < 1.25 pc']
  
  multiplot                     ;,/doyaxis
  ;; Third Panel *****************************
  
  ;; Compute PL fit
  mass = 10.^phys[i3].mass[0]
  mmax = m_max
  alpha = OMNI_MFN_PL(mass, mmax, mmin)
  
  OMNI_PLOT_CUTS, mass, /last
  
  al_legend,/top,/left,box=0,textcolor='red',$
            [string(cgSymbol('alpha'),alpha,format="(A0,' = ',F0.2)"),$
             'Cloud','R > 1.25 pc']
  
  
  myps,/done,/mp

  
  ;;******************************************************************
  ;; Number Density Cuts
  print,''
  print,'--------------------------------'
  print,'Number Density'
  
  ires = where(finite(s.rad),nres)
  iunr = where(~finite(s.rad),nunr)
  nden = phys[ires].nden[0]
  
  myps,'./masses/plots/astrocut_nden.eps',xsize=10
  print,!x.omargin,!y.omargin
  !x.omargin = [-7,1]
  !y.omargin = [-1,-3]
  multiplot,[2,1],xgap=0.002
  sunsym = cgSymbol('sun')
  
  i1 = where(nden LE alog10(750), n1)
  i2 = where(nden GT alog10(750) AND nden LE alog10(1.d4), n2)
  i3 = where(nden GT alog10(1.d4), n3)
  nt = n1+n2+n3+nunr
  
  
    ;; ;; Second Panel *****************************
  
  ;; Compute PL fit
  mass = 10.^phys[i1].mass[0]
  mmax = m_max
  alpha = OMNI_MFN_PL(mass, mmax, mmin)
  
  OMNI_PLOT_CUTS, mass, /first
  
  al_legend,/top,/left,box=0,textcolor='red',$
            [string(cgSymbol('alpha'),alpha,format="(A0,' = ',F0.2)"),$
             'Cloud','n < 750 cm!u-3!n']
  
  multiplot                     ;,/doyaxis
  ;; Third Panel *****************************
  
  ;; Compute PL fit
  mass = 10.^phys[i2].mass[0]
  mmax = m_max
  alpha = OMNI_MFN_PL(mass, mmax, mmin)
  
  OMNI_PLOT_CUTS, mass, /last
  
  al_legend,/top,/left,box=0,textcolor='red',$
            [string(cgSymbol('alpha'),alpha,format="(A0,' = ',F0.2)"),$
             'Clump','750 cm!u-3!n < n < 10!u4!n pc']
  
  
  myps,/done,/mp
  
  
  ;;******************************************************************
  ;; Linewidth Cuts
  print,''
  print,'--------------------------------'
  print,'Linewidth'
  
  idg = where(v.vlsr GE -500 AND ~v.inflw, ndg) ; Has DG tracer & finite LW
  ico = where(v.vlsr LT -500 AND ~v.inflw, nco) ; Has 13CO VLSR & finite LW
  inf = where(v.inflw, ninf)                    ; Infinite LW
  cn = [v[idg].cnum,v[ico].cnum] ; Order of CNUM for concatenation
  lw = [v[idg].lw,v[ico].grs.lw] ; Concatenated linewidths
  si = sort(cn)                  ; Indices of lw sorted by CNUM for use 
  lw = lw[si]                    ; Linewidths back in proper order
  
  myps,'./masses/plots/astrocut_linw.eps',xsize=10
  print,!x.omargin,!y.omargin
  !x.omargin = [-7,1]
  !y.omargin = [-1,-3]
  multiplot,[3,1],xgap=0.002
  sunsym = cgSymbol('sun')
  
  i1 = where(lw LE 1.5, n1)
  i2 = where(lw GT 1.5 AND lw LE 3.0, n2)
  i3 = where(lw GT 3.0, n3)
  nt = n1+n2+n3+ninf
  print,nt  
  
  ;; Compute PL fit
  mass = 10.^phys[i1].mass[0]
  mmax = m_max
  alpha = OMNI_MFN_PL(mass, mmax, mmin)
  
  OMNI_PLOT_CUTS, mass, /first
  
  al_legend,/top,/left,box=0,textcolor='red',$
            [string(cgSymbol('alpha'),alpha,format="(A0,' = ',F0.2)"),$
             'Quiet','lw < 1.5 km s!u-1!n']
  
  multiplot                     ;,/doyaxis
  ;; Second Panel *****************************
  
  ;; Compute PL fit
  mass = 10.^phys[i2].mass[0]
  mmax = m_max
  alpha = OMNI_MFN_PL(mass, mmax, mmin)
  
  OMNI_PLOT_CUTS, mass
  
  al_legend,/top,/left,box=0,textcolor='red',$
            [string(cgSymbol('alpha'),alpha,format="(A0,' = ',F0.2)"),$
             'Medium','1.5 km s!u-1!n < lw < 3.0 km s!u-1!n']
  
  multiplot                     ;,/doyaxis
  ;; Third Panel *****************************
  
  ;; Compute PL fit
  mass = 10.^phys[i3].mass[0]
  mmax = m_max
  alpha = OMNI_MFN_PL(mass, mmax, mmin)
  
  OMNI_PLOT_CUTS, mass, /last
  
  al_legend,/top,/left,box=0,textcolor='red',$
            [string(cgSymbol('alpha'),alpha,format="(A0,' = ',F0.2)"),$
             'Loud','lw > 3.0 km s!u-1!n']
  
  
  myps,/done,/mp
  
  
  ;;******************************************************************
  ;; Visual Extinction Cuts
  print,''
  print,'--------------------------------'
  print,'Visual Extinction'
  
  s[iunr].rad = min(s[ires].rad)
  
  Av = 42.4 * s.(iflux) / ( s.rad * s.rad / $
                            14. / 14. )
  Av = Av[ires]
  
  myps,'./masses/plots/astrocut_avis.eps',xsize=10
  print,!x.omargin,!y.omargin
  !x.omargin = [-7,1]
  !y.omargin = [-1,-3]
  multiplot,[2,1],xgap=0.002
  sunsym = cgSymbol('sun')
  
  i1 = where(Av LE 8., n1)
  i2 = where(Av GT 8., n2)
  nt = n1+n2+nunr
  print,nt
  
  ;; Second Panel *****************************
  
  ;; Compute PL fit
  mass = 10.^phys[i1].mass[0]
  mmax = m_max
  alpha = OMNI_MFN_PL(mass, mmax, mmin)
  
  OMNI_PLOT_CUTS, mass, /first
  
  al_legend,/top,/left,box=0,textcolor='red',$
            [string(cgSymbol('alpha'),alpha,format="(A0,' = ',F0.2)"),$
             'Low Column','A!dV!n < 8 mag']
  
  multiplot                     ;,/doyaxis
  ;; Third Panel *****************************
  
  ;; Compute PL fit
  mass = 10.^phys[i2].mass[0]
  mmax = m_max
  alpha = OMNI_MFN_PL(mass, mmax, mmin)
  
  OMNI_PLOT_CUTS, mass, /last
  
  al_legend,/top,/left,box=0,textcolor='red',$
            [string(cgSymbol('alpha'),alpha,format="(A0,' = ',F0.2)"),$
             'High Column','A!dV!n > 8 mag']
  
  
  myps,/done,/mp
  
  
  
  ;;******************************************************************
  ;; Galactocentric Azimuth Cuts
  print,''
  print,'--------------------------------'
  print,'Galactocentric Azimuth'
  
  phigal = atan(-xgal, -ygal) * !const.rtod
  
  myps,'./masses/plots/astrocut_gcaz.eps'
  
  cgPlot,phigal,pvec.stat.duse[0]/1.d3,psym=16,symsize=0.4,$
         xtit='Galactocentric Azimuth  [degr]',$
         ytit='Heliocentric Distance  [kpc]',charsize=1.0
  
  vline,50,color='red'
  vline,/h,6.5,color='red'
  
  linfit_plot,/robust,/legend,/top,/left,color='cyan',$
              phigal,pvec.stat.duse[0]/1.d3
  
  rho = correlate(phigal,pvec.stat.duse[0]/1.d3)
  al_legend,/bottom,/right,box=0,[cgSymbol('rho')+' = '+$
                                  string(rho,format="(F0.3)")]
  
  myps,/done
  
  ;;******************************************************************
  ;; Galactic Longitude Cuts
  print,''
  print,'--------------------------------'
  print,'Galactic Longitude'
  
  myps,'./masses/plots/astrocut_glon.eps',xsize=10
  print,!x.omargin,!y.omargin
  !x.omargin = [-7,1]
  !y.omargin = [-1,-3]
  multiplot,[2,1],xgap=0.002
  sunsym = cgSymbol('sun')
  
  i1 = where(pvec.glon LE 90., n1)
  i2 = where(pvec.glon GT 90., n2)
  nt = n1+n2
  print,nt
  
  ;; Second Panel *****************************
  
  ;; Compute PL fit
  mass = 10.^phys[i1].mass[0]
  mmax = m_max
  alpha = OMNI_MFN_PL(mass, mmax, mmin)
  
  OMNI_PLOT_CUTS, mass, /first
  
  al_legend,/top,/left,box=0,textcolor='red',$
            [string(cgSymbol('alpha'),alpha,format="(A0,' = ',F0.2)"),$
             'Quadrant I','l < 90 degr']
  
  multiplot                     ;,/doyaxis
  ;; Third Panel *****************************
  
  ;; Compute PL fit
  mass = 10.^phys[i2].mass[0]
  mmax = m_max
  alpha = OMNI_MFN_PL(mass, mmax, mmin)
  
  OMNI_PLOT_CUTS, mass, /last
  
  al_legend,/top,/left,box=0,textcolor='red',$
            [string(cgSymbol('alpha'),alpha,format="(A0,' = ',F0.2)"),$
             'Quadrants II & III','l > 90 degr']
  
  myps,/done,/mp
  
  ;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  ;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  ;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  ;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  ;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  ;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  ;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  ;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  ;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  ;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  ;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  ;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  ;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  ;; Footnote junk to say something about the state of the world.
  time = systime(1) - t_start
  hour = time/3600.
  min  = (time - floor(hour)*3600)/60.
  sec  = time - floor(hour)*3600 - floor(min)*60
  message,'This run had sampletd = '+$
          string(sampletd,format="(I0)")+'.',/inf
  message,'For '+string(nreal,format="(I0)")+' realizations of the mass '+$
          'function, the routine took '+string(floor(hour),format="(I0)")+$
          ' hours, '+string(floor(min),format="(I2)")+' minutes, and '+$
          string(round(sec),format="(I2)")+' seconds.',/inf
  RETURN                        ; of the Jedi
END
